Transfer Function Identification using 
Orthogonal Fourier Transform Modeling Functions 

Eugene A. Morelli 1 

NASA Langley Research Center, Hampton, Virginia, 23681 


A method for transfer function identification, including both model structure 
determination and parameter estimation, was developed and demonstrated. The approach 
uses orthogonal modeling functions generated from frequency domain data obtained by 
Fourier transformation of time series data. The method was applied to simulation data to 
identify continuous-time transfer function models and unsteady aerodynamic models. 
Model fit error, estimated model parameters, and the associated uncertainties were used to 
show the effectiveness of the method for identifying accurate transfer function models from 
noisy data. 




Nomenclature 

Cl 

= 

nondimensional aerodynamic lift coefficient 


= 

expectation operator 


= 

Fourier transform operator 

Im[ • ] 

= 

imaginary part 

j 

= 

imaginary number = «s/-F 

J 

= 

cost function 

N 

= 

number of data points 

PSE 

= 

predicted squared error 

**[■] 

= 

real part 

s 

= 

Laplace variable 

t 

= 

time, s 

T 

= 

maneuver length, s 

a 

= 

angle of attack, rad or deg 

e 

= 

parameter vector 

i 

= 

covariance matrix 

CO 

= 

frequency, rad/s 

superscripts 



T 

= 

transpose 

* 

= 

estimate 

~ 

= 

transform 

(0 

= 

i th time derivative 

-i 

= 

matrix inverse 
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subscripts 

R 

C 

S 

U 


real 

complex 

quasi-steady 

unsteady 


I. Introduction 

JN the excellent textbook entitled Numerical Methods That (Usually) Work , author Forman Acton writes 1 : 

A favorite form of lunacy among aeronautical engineers produces countless attempts to decide what 
differential equation governs the motion of some physical object, such as a helicopter rotor. In the 
archetypical experiment, one measures both the input to the blade and its output - that is, motion. 
Arguing that the system is a (( black box ” to which they have the input and the output, they now wish 
to compute the transfer function, another name for the differential equation. The difficulty here is that 
many quite different differential equations can give rise to nearly the same outputs from identical 
inputs. Thus discriminating the (( true ” equation from its brethren is next to impossible. 

In general, a variety of mathematical model structures can be used to fit a particular set of input-output data 
within a selected error tolerance. For the transfer function model form, this translates to the well-known fact that 
high-quality model fits to measured data can be obtained using many different orders for the numerator and 
denominator. The critical difference among such models, not mentioned in the quote above, is that higher-order 
models do not predict as well for new data, because excess model structure generally fits noise, which is not 
repeatable. On the other hand, a model with insufficient model structure will not capture important characteristics 
exhibited in the data, and this also adversely affects prediction accuracy. 

In much of the literature on transfer function modeling 2 * ' 6 , the model structure (equivalently, the order of the 
transfer function numerator and denominator) is assumed to be known, so that the modeling problem is reduced to 
parameter estimation. Extensive work has been done on the so-called “order determination problem” or “variable 
selection problem” 7 ' 8 , which is equivalent to determining the structure of the transfer function model, mainly for 
discrete-time models, but also for continuous-time models 2,8 . The quote from Professor Acton’s book highlights the 
need for more work in this area. Transfer function models are important in aerodynamic modeling, aircraft 
simulation, control system design, actuator modeling, and many other applications, so it is important that models 
identified from measured data be accurate, with good prediction capability. The model must also be identifiable 
based on a defined criterion that distinguishes various model forms that are essentially equivalent in terms of model 
fit to the data. 

In this paper, a method is developed for identifying transfer function model structure (numerator and 
denominator orders) as well as the associated unknown model parameters (transfer function coefficients) and 
associated uncertainties, based on measured data alone. The method uses orthogonal modeling functions derived 
from finite Fourier transforms of the measured input-output data. The development here is focused on single-input 
single-output (SISO) transfer function identification, although the approach could also be applied to multi-input 
multi-output (MIMO) transfer function matrix identification. 

All of the experiment design, data analysis, and modeling tasks included in this work were done using system 
identification software written in MATLAB®, called System IDentification Programs for AirCraft, or SIDPAC 9 . 
SIDPAC is bundled with Ref. [9], and is therefore publicly available. 

The next section describes the method used to identify transfer function model structure and estimate the 
associated model parameters and uncertainties. This section also includes discussion of a statistical modeling metric 
used for model structure determination, along with a modulating function approach that can be used to handle 
general non-zero endpoint conditions, and a discussion of real-time implementation. In the Examples section, the 
modeling approach is applied to a simulation problem where the true model is known and measurement noise is 
added to the output. Following that, the modeling technique is applied to identify an aerodynamic transfer function 
for unsteady aerodynamics. Finally, concluding remarks summarize and discuss the results. 
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II. Method 

The method for both model structure determination and parameter estimation is based on first transforming 
sampled time-domain data into the frequency domain. This transformation is done using a high-accuracy 
implementation of the finite Fourier transform at arbitrarily -chosen frequencies 9,10 . Next, the complex-valued 
frequency domain data is orthogonalized and processed using a statistical modeling metric for efficient and 
automated model structure determination. Finally, the selected orthogonal model terms are transformed back into 
physically-relevant model terms, and estimated model parameter values and uncertainties are calculated for the final 
model. 

A. Finite Fourier Transform 

The finite Fourier transform is the analytical tool used to transform measured data from the time domain into the 
frequency domain. For an arbitrary scalar signal x(V) on the time interval [0, T] , the finite Fourier transform is 
defined by 


J r [v(/ L )] = x(co) = J x(t^)e j wt dt 


( 1 ) 


which can be approximated by 


xWvAt^xfte-W' ( 2 ) 

i = o 

where x(i) = x(iAt ) and At is a constant sampling interval. The summation in Eq. (2) is defined as the discrete 
Fourier transform, 


X M = Z *(0 e ~ jmiAt ( 3 ) 

i = 0 

so that the finite Fourier transform approximation in Eq. (2) can be written as 

x(co)nX(co) At (4) 

Eq. (4) represents a simple Euler approximation to the finite Fourier transform in Eq. (1). Cubic interpolation can be 
applied to the sampled data x(i) to improve the numerical accuracy of the transformation, and arbitrary frequencies 

can be selected for the finite Fourier transform calculation. With these modifications, the finite Fourier transform in 
Eq. (1) can be computed numerically for arbitrary frequencies, with accuracy on the order of the numerical precision 
of the computer. The calculations are explained in Refs. [9] and [10], and were used to compute the finite Fourier 
transform in Eq. (1). 

When transforming time-domain data into the frequency domain, bias and drift are always removed from each 
measured time series prior to applying the Fourier transform. This is to avoid spillage from large low-frequency 
components which can pollute the frequency -domain data at low frequencies of interest 9 . 

When a time series has non-zero endpoints, the Fourier transform of the time derivative of the time series 
requires endpoint correction terms. Using integration by parts, the finite Fourier transform of the time derivative of 
x(V) is 


= x(T)e i <oT - v(0) + jco v( of) 


( 5 ) 


Similarly, 
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( 6 ) 


f[x(t)~\ = x(T) e~ ja>T - x(0) + jcoT\_x{t)\ 

= x(T)e~ jaT -i(0) + jco[x(T)e- J(oT -Jt(0) ]+ (jcof x(co) 

and so on for higher-order derivatives. When the magnitude of the data and time derivatives at the endpoints are 
small, e.g., for a flight maneuver that starts and ends at the same steady trim condition, then this effect is small, and 
can be ignored. However, the effect cannot be ignored in general, and will involve correction terms requiring up to 

{n — l) th order derivatives at the endpoints for the Fourier transform of n th order derivatives of the time series. 

B. Orthogonal Fourier Transform Modeling Functions 

In past work 911 ' 13 , real data vectors were orthogonalized to effectively and accurately determine model structure 
and associated model parameters for general nonlinear aerodynamic models. For transfer function modeling in the 
frequency domain, the data are complex. This requires slightly different data handling to achieve the 
orthogonalization necessary for model structure determination, but otherwise, the same approach used for real data 
can be applied. 

The necessary modification for complex data is very simple, and is related to what is called the Euclidean angle 
between two complex vectors 14 . For an arbitrary complex vector x c with N generally complex elements, a real 
isometric vector x R with 2N elements can be assembled by stacking the real and imaginary parts of the complex 
elements of x c : 


M* c ] 

Jm[x c ] 


(7) 


Treating the real vector x R as the new data vector, all of the orthogonalization techniques and model structure 
determination techniques developed previously for real data can be applied without modification. The penalty 
associated with applying the techniques to complex data is simply a doubling of the length of the data vectors 
involved. The next two subsections give a brief overview of the orthogonal function modeling technique. Further 
details can be found in Refs. [9], [1 1 ]-[13], which are the original sources for the technique. 

7. Generating Orthogonal Modeling Functions 

Mutually orthogonal vectors P\> P 2 > •••> Pn Q can generated from real vectors g lf g 2 , • • €n 0 us i n g an 

orthogonalization procedure analogous to Gram-Schmidt orthogonalization for vectors. The process begins by 
choosing one of the f • as the first orthogonal vector, e.g., p l =^ l . In general, any f • can be chosen as the first 

orthogonal vector, without any change in the procedure. The orthogonal vectors are then generated sequentially as 


j ■ 

= $j-Y.rkjPk j = 2,3,...,n a 


k = 1 


( 8 ) 


The y k j for k — 1, 2, — l are scalars determined by multiplying both sides of Eq. (8) by pi , then invoking the 

mutual orthogonality of the p k , k — 1,2, ..., j , and solving for y k j 


y kj = 4^- j = 2,3,...,n 0 k =1,2,..., j -1 (9) 

PkPk 

It can be seen from Eqs. (8)-(9) that each orthogonal vector can be expressed exactly in terms of a linear 
expansion of the original vectors. The orthogonal vectors are generated sequentially by orthogonalizing the original 
vectors with respect to the orthogonal vectors already computed. Each orthogonal vector can be considered an 
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orthogonalized version of an original vector. Because the original vectors g lf g 2 , • • € Ho are g enera l functions of 

measured data, these vectors can be called original functions, and the corresponding orthogonal vectors 
Pi y p 2 , • • p n can be called orthogonal functions. 

2. Conversions Between Orthogonal Functions and Original Functions 

If the pj vectors and the gj vectors are arranged as columns of matrices P and X , respectively, and the y k j 

are elements in the k th row and / h column of an upper triangular matrix G with ones on the diagonal, 


then 


and 


p-[pi pi - 


*=[f> fe - ij 


1 

0 


y 12 

i 

o 


r is 

y 23 

i 


rm o 
y 2 n a 

y^ a 


0 0 0 ••• 1 


(10a) 


(10b) 


X = PG 


(ID 


P = XG~ l 


( 12 ) 


The columns of G 1 contain the coefficients for expansion of each orthogonal function (column of P ) in terms 
of an exact linear expansion in the original functions (columns of X ). Equation (12) can be used to express each 
orthogonal function in terms of the original functions. The manner in which the orthogonal functions are generated 
allows them to be decomposed without ambiguity into an expansion of the original functions, which have physical 
meaning. 

The next section describes orthogonal function modeling in the frequency domain for the transfer function 
identification problem. 

C. Transfer Function Model Structure Determination and Parameter Estimation in the Frequency Domain 

The transfer function model for a response y(f) to an input u(t ) for zero initial conditions can be written as: 

y( s ) _ c q + c 1 s + c 2 s 2 +... + c m s m 
u(s) 1 + d x s + d 2 s 2 + . .. + d n s n 

where the c 0 ,c 1 ,c 2 ,...,c m and d } ,d 2 ,. ..,d n coefficients are model parameters, and the selection of terms for the 
numerator and denominator (which also determines the integer values of m and n) defines the model structure. 
Setting s = jco gives the equivalent model in terms of the Fourier transform, 


y[co) _ c 0 +c 1 ( 7 ^) + c 2 ( 7 ft>) 2 +... + c m ( jco^ 
11 ( (U ) l + d l [jco) + d 2 (jco) + ... + d n (jco) 


Re-arranging, 
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y (co) = c 0 u (co) + q (jco) u (co) + c 2 ( 7<w) 2 a (®) + . . . + c m (;&>/” w (<w) 
- d \ (ja)y(co) - d 2 (jco) 2 y(co)-...-d n (jco) n y(co) 


( 15 ) 


Allowing for noise s on the measured output z, 
z(co) = y(a>) + £ 

z(a>) = c 0 u(co) + c 1 (jeo)u(eo) + c 2 (j&>) 2 u(cd) + ... + c m ( jco) m u(co) (16) 

- d \ ( jco) z ( co) - d 2 ( jcof z (co) — ... — d n ( jcof z(co) + e' 

For identification purposes, this is an equation-error linear regression in the frequency domain with complex- 
valued regressors. If the complex regressors for selected frequencies are converted to real vectors using Eq. (7), and 
orthogonalized as described in the previous section, then the form of the orthogonal function model is 

z = c/| p l +a 2 p 2 + ... + a n p n +e (17) 


where z is a 2M-dimensional vector for the measured response in the frequency domain, 

Re\z(co)\ 

z— L J 


(18) 


and the quantity z (co) is the complex vector of Fourier transforms for M selected frequencies in the frequency 
vector co . The response z is modeled in terms of a linear combination of n s selected orthogonal functions 
Pj, j = 1,2 where n s < n Q , n Q =n + m + 1 , which allows for the possibility that not all of the model terms in 
Eq. (16) may be needed in the final model. Each pj is an 2M-dimensional vector which depends on input or output 
data and the analysis frequencies. For example, 


p i=fi = 


Re^u(co)~) 

Im^u(a))~) 


(19) 


The cij , j = 1,2 ,...,n s in Eq. (17) are constant model parameters to be determined, and £ denotes the modeling error 
vector. 

Now define an M xn s matrix P , 


and let 


P = [Pl’P2’->Pn s ] 


a = ^c 




( 20 ) 


( 21 ) 


Equation (17) can then be written as a standard least squares regression problem, 

z = Pa + e (22) 
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The error vector s is to be minimized in a least squares sense. The goal is to determine a that minimizes the least 
squares cost function 

J =^e T e = ^(z-Paf (z-Pa) (23) 

The parameter vector estimate a that minimizes this cost function is computed as 9 

a = [ pT p]~ lpT z (24) 

The estimated parameter covariance matrix is 9 


*a=E 




= cr 2 (p T P 




(25) 


9 

where the fit error variance cr can be estimated from the residuals 

v =z- Pd 


(26) 


using 


c^ 2 = 


(2 M-n s ) 


(z-Paf (z-Pd) 


T 

V V 


(2 M-n s ) 


(27) 


Parameter standard errors are computed as the square root of the diagonal elements of the 2^ matrix from Eq. (25), 
using cr from Eq. (27). The identified model output y is 

y = Pa (28) 


When the modeling functions are mutually orthogonal, an appropriate model structure can be determined in a 
straightforward manner, because the explanatory capability of each orthogonal modeling function is completely 
distinct from all of the others. This decouples the least squares modeling problem, as will be shown now. 

For mutually orthogonal modeling functions, 


pjpj = 0 , i*j , i,j ~ 1, 2, .... n () 


(29) 


so that for P = |~pj p 2 ... /?„ J, P r P is a diagonal matrix with the inner product of the orthogonal functions 

on the main diagonal. Using Eqs. (20)-(21) and (29) in Eq. (24), the / h element of the estimated parameter vector a 
is given by 


d j={p T j z )l(p T i Pj) j = \X...,n Q 


Using Eqs. (20)-(21), (29), and (30) in Eq. (23), 


y = i 
2 


zTz -2L(p]z$ l(p T jPj] 

7=1 


(30) 


( 31 ) 
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Equation (31) shows that when the modeling functions are orthogonal, the reduction in the least-squares cost 
function resulting from including the term aj pj in the model depends only on the response variable data z and the 

added orthogonal modeling function pj . The least-squares modeling problem is therefore decoupled, which means 

each orthogonal modeling function can be evaluated independently in terms of its ability to reduce the least-squares 
model fit to the data, regardless of which other orthogonal modeling functions are already selected for the model. 
When the modeling functions pj are instead a non-orthogonal function set, the least-squares problem is coupled, 

and iterative analysis is required to find the subset of modeling functions for an adequate model structure. 

The orthogonal modeling functions to be included in the model are chosen to minimize predicted squared error, 
PSE, defined by 15 


PSE = 


_ ( z-Pa) T ( z-Pa ) 


2 M 


+ cr z 


2 M 


(32) 


or 


pSE =—+^Lx— 

M 2 M 


(33) 


where n s is the number of orthogonal functions selected for the model, n s <n Q . 

The constant cj^ nax is the upper-bound estimate of the squared error between future data and the model, i.e., the 
upper-bound mean squared error for prediction cases. The upper bound is used in the model over-fit penalty term to 
account for the fact that PSE is calculated during the model structure determination stage, when the model structure 

is not correct. Using <J^ nax in the model complexity term also results in the PSE metric conservatively estimating 
the squared error for prediction cases. 

A simple estimate of <J^ nax that is independent of the model structure can be obtained by considering cjf nax to be 
the residual variance estimate for a constant model equal to the mean of the measured response values, 


2 


1 


2 M 


max oa/r t Z ] 


2M -1 


i = 1 


, 2 M 


i=\ 


(34) 


The PSE in Eq. (33) depends on the mean squared 
error, j/m , and a term proportional to the number of 
terms in the model, n s . The latter term prevents over- 
fitting the data with too many model terms, which is 
detrimental to model prediction accuracy 9,15 . While the 

mean squared error J/M must decrease with the 
addition of each orthogonal modeling function to the 
model (by Eq. (31)), the over-fit penalty term 

<J max n s/{2M) must increase with each added model 
term ( n s increases). Introducing the orthogonal 
modeling functions into the model in order of most 
effective to least effective in reducing the mean squared 

error (quantified by Pj ) for the / h 


x id 5 



Orthogonal Function Number 


Figure 1. Model structure determination using 
orthogonal functions and PSE 


orthogonal modeling function) results in the PSE metric 
always having a single global minimum. 

Figure 1 depicts this graphically, using actual modeling results. The figure shows that after the first 6 modeling 
functions, the added model complexity associated with the next most effective orthogonal modeling function is not 
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justified by the associated reduction in mean squared error. This point is marked by minimum PSE , which defines 
an adequate model structure with good predictive capability. Ref. [15] contains further statistical arguments and 
analysis for the form of PSE given in Eq. (33), including justification for its use in modeling problems. 

Using orthogonal functions to model the response variable makes it possible to evaluate the merit of including 
each modeling function individually , using the predicted squared error PSE. The goal is to select a model structure 
with minimum PSE, and the PSE always has a single global minimum for orthogonal modeling functions ranked 
according to their effectiveness in reducing the mean squared error. This makes the model structure determination a 
well-defined and straightforward process that can be (and was) automated. For transfer function model structure 
determination, selecting appropriate orthogonal functions is equivalent to choosing the appropriate numerator and 
denominator terms in the model. 

After the model structure is determined by selecting orthogonal modeling functions for minimum PSE, the 
identified model output is 


y = Pd (35) 

where the P matrix now includes only the n s orthogonal functions selected in the model structure determination, 
n s <n 0 . In general, the retained orthogonal functions will not necessarily be consecutive or the first in the 
orthogonalization sequence of Eq. (8). Each retained orthogonal function can be decomposed without error into an 

expansion of the original functions, using the columns of G 1 in Eq. (12) corresponding to the retained orthogonal 
functions. Common terms are combined using double precision arithmetic to arrive finally at a model using only 
original functions. Terms that contribute less than 0.1 percent of the final model root-mean- square magnitude are 
dropped. 

From Eq. (8), it is clear that the orthogonalization is dependent on the order in which the original functions are 
processed. If a selected orthogonal function happens to be late in the order, the result is that final conversion from 
the selected orthogonal function back to original functions will involve small parts of original functions not needed 
in the final model. To avoid this, orthogonalization is done initially to determine the order of most effective to least 
effective orthogonal functions, then the associated original functions are re-ordered to match this order of 
importance, and the orthogonalization is repeated. This simple process results in the minimum number of model 
terms for the final model in terms of original functions. 

The final form of the model is a sum of ordinary functions, with associated model parameter estimates and 
uncertainty. For the transfer function modeling problem, the model consists of selected terms from the complete set 
of postulated terms in Eq. (16). 

D. Using Modulating Functions to Avoid Endpoint Corrections 

The finite Fourier transform of time derivatives involves endpoint corrections when the time series endpoints 
and/or derivatives are non-zero, as shown earlier in Eqs. (5) and (6). These endpoint corrections can be avoided 
using modulating functions, as first described by Shinbrot 16 . The general idea is to multiply the differential equation 
associated with Eq. (13) by a suitably-chosen modulating function, then use integration by parts to transfer the 
derivatives from the data to the modulating function, incurring endpoint terms that involve the modulating functions 
and its time derivatives evaluated at the endpoints. A smooth differentiable function (j){t ) with the value of the 

function and its first n — 1 time derivatives equal to zero at the endpoints is called an n th -order modulating function, 

/ k \o) = / k \T) = 0 k = 0,l,...,n-l (36) 


where (jP* (t) denotes the first time derivative of </;(?) , and so on. 

Pearson 17 developed a Fourier-based modulating function. The Fourier-based modulating function for frequency 
co is 


4>(t) = e~ ja>t (l-e~ jm A n co a = 2 n/T 


( 37 ) 
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The function ^ i>(t ) in Eq. (37) satisfies the requirements of Eq. (36) for an n th -order modulating function. This 
modulating function can also be written as 


Mt) = t b k e- j ^ ka °> b k = (-if t n! - 


(38) 


which shows that the modulating function (p(t) is a summation of complex sinusoidal functions of the same type 
used in the finite Fourier transform of Eq. (1). 

Starting with the differential equation associated with Eq. (16), 


Z = CqU + C^U W + c 2 u + . . . + c m u^ - d x z^ - d 2 z - ... - d n z^ + s’ 
Multiplying the equation by (p(t) , and integrating over the time period [0, T ] , 


(39) 


(pzdt-c oj^ (pudt + c x (pu^ dt + c 2 ^ (pu^ dt + + (pu^dt 

</>z^dt — d 2 ^ <pz^ dt-...-d n (pz^dt + s' 


where s' represents the complex- valued equation error for the modulated equation. Using integration by parts, and 
invoking the properties of an n ih -order modulating function from Eq. (37), 


(pz dt - c 0 J o (pudt-c x <p( 1 ^ udt + c 2 ^ (p( 2 ^ udt - ... 

+ ^1^ <p(^ zdt-d 2 j^ (p^ zdt + ... + d n (~l) n+1 (p^ zdt + s' 


(40) 


Each term in Eq. (40) is a complex-valued modulation of the input-output data for frequency co . Carrying out 
this process for a vector of frequencies co provides frequency-domain data of the same type obtained earlier using 
the finite Fourier transform. The difference is that using the finite Fourier transform required assuming that the 
endpoint conditions were close to zero, so that the endpoint terms in Eqs. (5) and (6) would be negligible. Using the 
modulating function approach, that assumption is not required, because the modulation process removes the 
non-zero endpoint terms analytically. The modulating function approach therefore allows the use of input-output 
time series data with arbitrary endpoint conditions. 

At this point, the model structure determination and parameter estimation proceed in exactly the same manner as 
described earlier in Sections II.B. and II. C. Because the modulating function in Eq. (38) uses complex exponential 
functions, the required modulation integrals in Eq. (40) can be carried out using the same high -accuracy, arbitrary- 
frequency finite Fourier transform calculation used earlier 9,10 . 

E. Real-Time Implementation 

Transfer function identification in the frequency domain can also be implemented in real time, if the calculations 
described earlier are applied to frequency-domain data at selected times, typically at regular time intervals. Because 
the method uses frequency-domain data from the Fourier transform, it is necessary to have values of the Fourier 
transform available at any time. Fourier transforms can be computed in real time using a recursive Fourier 
transform, as described next. 

For a given frequency co , the discrete Fourier transform in Eq. (3) at time i At , denoted by X t (co) , is related to 
the discrete Fourier transform at time (i - 1 )At by 


10 

American Institute of Aeronautics and Astronautics 



(41) 


X i (co) = X i _ l (co) + x(i)e- ja > iAt 


where 


e ~j(oiAt _ e ~jcoAt e ~jco(i-\)At 


(42) 


The quantity e ^ At is constant for a given frequency co and constant sampling interval At . It follows that the 
discrete Fourier transform can be computed for a given frequency at each time step using one addition in Eq. (41) 

and two multiplications - one in Eq. (42) using the stored constant e~^ At for frequency co , and one in Eq. (41). 
There is no need to store the time-domain data in memory when computing the discrete Fourier transform in this 
way, because the data for each sample time is processed immediately. Time-domain data from the past can be used 
in all subsequent analysis by simply continuing the recursive calculation of the Fourier transform. In this sense, the 
recursive Fourier transform implements a memory for spectral information in the data. More data from more 
maneuvers improves the quality of the data in the frequency domain without increasing memory requirements to 
store it. Furthermore, the Fourier transform is available at any time i At . The approximation to the finite Fourier 
transform is completed using Eq. (4). 

The recursive computation of the Fourier transform does not use a Fast Fourier Transform (Err) algorithm 18 , 
but rather computes the Fourier transform in a straightforward way, using Eqs. (41) and (42). The recursive Fourier 
transform would therefore be comparatively slow, if the entire frequency band up to the Nyquist frequency l/(2 At) 

were of interest. However, the dynamics for many physical systems (such as rigid-body aircraft dynamics, rotor 
dynamics, or actuator dynamics) lie in a rather narrow frequency band of approximately 0.01-2.0 Hz. Because the 
frequency band is limited, it is efficient to compute the discrete Fourier transform using Eqs. (41) and (42) (which 
represents a recursive formulation of Eq. (3)) for only selected frequencies. With this approach, it is possible to 
select closely-spaced fixed frequencies for the Fourier transform and the subsequent modeling, and still do the 
calculations efficiently. 

Using a limited frequency band for the Fourier transformation confines the data analysis to the frequency band 
where the dynamics lie, and automatically filters wide band measurement noise or other dynamics (such as structural 
responses) outside the frequency band of interest. These automatic filtering features are important for real-time 
applications, where instrumentation error corrections and noise filtering would require additional computational 
resources. 

In past work using the recursive Fourier transform for stability and control parameter estimation on fighter 
aircraft and subscale aircraft 9,19,20 , frequency spacing of 0.04 Hz on an interval of approximately [0.1-2] Hz was 
found to be adequate. Finer frequency spacing requires slightly more computation, but was found to have little 
impact on the results. When the frequency spacing is very coarse, there is a danger of omitting important frequency 
components, and this can lead to inaccurate parameter estimates. In general, a good rule of thumb is to use 
frequencies evenly spaced at 0.04 Hz over the bandwidth for the dynamic system. For good results, the bandwidth 
should be limited to the frequency range where the signal components in the frequency domain are at least twice the 
amplitude of the wide band noise level. However, the modeling is robust to these design choices, so the selections 
to be made are not difficult. 

The recursive Fourier transform update need not be done for every sampled time point. Systematically skipping 
time points effectively lowers the sampling rate of the data prior to Fourier transformation. This saves computation, 
and does not have a significant adverse impact on the modeling results until the Fourier transform update rate gets 
below approximately 5 times the highest frequency of interest for the dynamic system. The model structure 
determination, parameter estimation, and covariance calculations described earlier can be done at any time, but 
would typically be done at 1 or 2 Hz, to save computations, except in cases of strong nonlinearity, damage, failure, 
or rapid changes in conditions. For these cases, the model identification update rate can be increased, at the cost of 
more computations. 

For real-time modeling, measured time series are high-pass filtered to remove the steady part of each signal. 
Then, each perturbation signal is transformed into the frequency domain using the recursive Fourier transform. The 
break frequency for the high-pass filter is set just below the lowest frequency selected for the modeling and Fourier 
transform. High -pass filtering can be implemented with a fourth-order Butterworth digital filter, for example. The 
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high-pass filtering prevents leakage from the relatively large spectral component at zero frequency, associated with 
the steady component of each signal, from polluting transformed data at low frequencies. 

F. Application to Frequency Responses 

Frequency responses are typically shown using the Bode plot, using separate plots of the magnitude in dB, and 
the phase angle in degrees, plotted against rotational frequency co using a semi-log scale. The frequency response 
data used to make the Bode plot consists of a vector of complex numbers associated with a vector of frequencies. 
Although it is simpler and more straightforward to use the finite Fourier transform or Fourier-based modulating 
functions, it is also possible to use frequency response data with the same modeling approach described in Sections 
II.B. and II. C. Only a slight modification is required because of the use of cross-spectral density and auto-spectral 
density in the computation of frequency responses. The analog of Eq. (16) is then 


S m H = c o Suu H + Cl (. jco ) S uu (co) + c 2 (jcof S uu (co) + ... +c m (jco) m S uu (co) 
~ d \ ( M) S zu (co)-d 2 ( jcof S zu (co)-...-d n (jcof S zu (co) + s' 


(43) 


where s' represents the complex- valued equation error. Using this form, all of the material on transfer function 
model structure determination and parameter estimation given earlier can be applied without modification. 

G. Application to State-Space Model Structure Determination in the Frequency Domain 

The model structure determination approach described in Sections II.B. and II. C. can also be applied to 
determine model structure for state-space models in the frequency domain. For this application, the procedure is 
identical to the approach taken in the time domain, described in Ref. [9], except that complex data in the frequency 
domain is first modified by the simple expression in Eq. [7] as a preliminary step. Beyond that, model structure 
determination procedures such as stepwise regression and orthogonal function modeling, which depend 
fundamentally on orthogonalization in function space, can be applied in exactly the same manner as in the time 
domain. This approach to state-space model structure determination in the frequency domain is more efficient and 
less subjective than the common practice 6 of iteratively estimating the parameters in candidate state-space model 
structures, then retaining or rejecting specific terms based estimated uncertainty measures or parameter 
insensitivities. 


III. Examples 

The transfer function identification procedure will be demonstrated using simulated data for known transfer 
functions. Application to identifying unsteady aerodynamic models will also be discussed and demonstrated. 

A. Transfer Function Identification from Simulated Data 

An optimized wide-band multisine excitation input 9,19-21 was applied to a known transfer function: 

y(d_ i+Q-5* (44) 

m(j) 1 + 0. 159s + 0.025s 2 

and the output was corrupted with 5 percent white Gaussian noise. The simulated time series data are shown in the 
top plot of Figure 2. The transfer function model in Eq. (44) could represent the short-period pitch rate dynamic 
response to elevator deflections for an airplane, but could also describe many other physical systems. 

The transfer function identification procedure was applied to the measured time-domain data, producing the 
modeling results given in Table 1. The results indicate accurate model structure determination and model parameter 
estimates, with appropriate and small error bounds for the estimated model parameters. Note that the algorithm 
identified the correct model structure (transfer function numerator and denominator orders) in addition to accurately 
estimating the model parameters. Maximum model order equal to 3 was used to generate the candidate orthogonal 
modeling functions in the frequency domain. 

The model fit to the time-domain output data is shown in the middle plot of Figure 2. The lower plot shows the 
model residuals in the time domain, which are the difference between the measured output and the model output. 
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The residuals are essentially white noise, indicating that the modeling process has successfully identified the 
deterministic functional dependencies in the data. 

Figure 3 shows the fit of the identified model to the measured output data in the frequency domain, along with 
the frequency-domain residual. The plots in Figure 3 use the modulus of the complex -valued frequency-domain 
quantities, which is the square root of the summed squares of the real and imaginary parts, 


y | = ^[^(y)] 2 + [/m(y)] 2 (45) 


The input multisine had component frequency content in the frequency band [0.3, 1.8] Hz, at evenly-spaced 
intervals of 0.3 Hz. The analysis frequencies were selected in the frequency band [0.1, 2] Hz at evenly-spaced 
intervals of 0.1 Hz. There is nothing critical about these selections, except that the analysis frequency band should 
encompass the frequency content of the multisine input. Many other choices work equally well. 

Figures 2 and 3, and the results in Table 1, show that the transfer function identification method produced an 
accurate transfer function model in both the time domain and the frequency domain, based on the input -output time 
series data alone. The modulating function approach works equally well, with results from applying that technique 
also given in Table 1. Model fit to the data using the modulating function approach was very similar to what is 
shown in Figs. 2 and 3. 
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Figure 2. Transfer function identification 
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Figure 3. Transfer function model fit in the frequency domain 
B. Unsteady Aerodynamic Modeling 

Unsteady aerodynamic models can be identified by applying the equation-error transfer function identification 
method using orthogonal Fourier transform modeling functions. For this application example, the analysis is 
confined to the response of lift coefficient C L to angle of attack a . The same procedure could be applied to other 
input-output pairs. 

A model for unsteady aerodynamic lift coefficient response to angle of attack can be built based on material in 
Refs. [22] and [23]. Using a differential equation to implement the unsteady aerodynamics, 


Q 

1 

II 

(46a) 

f] = —brj + a 

(46b) 

where C L ^ denotes the unsteady part of the lift coefficient and 7/ is an additional state associated with the unsteady 

aerodynamics. Model parameters a and b in general depend on the nominal angle of attack a Q . The differential 
equation model form above is equivalent to an indicial function model for unsteady lift, 

ill 

© 4. 

l 

<3- 

T 

S- 

'X 

(47a) 

~-arj--a j^e b ^ a (r) dr 

(47b) 


where the indicial function is postulated as the simple exponential function e bt . In transfer function form, the 
model for unsteady lift can be derived from Eqs. (46) as 
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as 


( 48 ) 



a s + b 


The unsteady aerodynamic model in Eqs. (46)-(48) describes an addition to the lift coefficient that depends on 
past values of a , but fades out with time. The b parameter quantifies duration of the unsteady effect, and the a 
parameter quantifies the magnitude of the unsteady effect. The unsteady aerodynamics C L ^ is added to the 

quasi-steady lift coefficient dependence on angle of attack, which is normally obtained from static wind tunnel tests 
or stability and control flight tests. 

For relatively small and slow variations in angle of attack, the transfer function for the lift coefficient due to 
angle of attack can be modeled as a constant C L , 


C L 

— = Cl u (49) 

a a 


The model in Eq. (49) implements the quasi-steady assumption that the lift coefficient depends only on the current 
angle of attack. The transfer function for the lift coefficient arising from arbitrary angle of attack variations is the 
sum of the quasi- steady and unsteady parts, 


Cl = C h , C K = ( C *» - fl ) 5 + c x a fr 

a a a (s + Z?) 

Converting to the form used for transfer function model identification, 

Cl = + ( Cz a ~ Q )^ _ c 0 + fig 

a 1 + b~ x s 1 + d\s 


(50) 


(51) 


This modeling represents the idea of identifying an aerodynamic transfer function for the influence of angle of 
attack on lift coefficient (cf. Ref. [24]), as opposed to the simple quasi-steady form shown in Eq. (49) that is 
typically derived from wind tunnel tests or stability and control flight tests. 

The top plot in Fig. 4 shows input-output data for simulated unsteady aerodynamics, using the model structure in 
Eq. (51). The angle of attack was varied using an optimized wide -band excitation input 9,19-21 and the simulated 
output was corrupted by 5 percent Gaussian noise. The lower plots in Fig. 4 and both plots in Fig. 5 show the model 
fit to the data in the time and frequency domains. The input multisine had component frequency content in the 
frequency band [0.04, 0.68] Hz, at evenly-spaced intervals of 0.04 Hz. The analysis frequencies were selected in the 
frequency band [0.04, 0.82] Hz at evenly-spaced intervals of 0.02 Hz. Again, there is nothing critical about these 
selections, except that the analysis frequency band should encompass the frequency content of the multisine input. 

Table 2 contains the transfer function identification results. As before, the transfer function identification 
procedure accurately identified both the model structure and the model parameters. Table 2 also contains the values 
of a , b , and C L used in the model of Eq. (51). 

This modeling approach is a practical implementation of the concept of generalizing stability and control 
derivatives to aerodynamic transfer functions, as described in theory by Etkin 24 , more than 55 years ago. Figure 6 
shows the identified aerodynamic transfer function for lift coefficient response to angle of attack from the last 
example. The plot indicates that the C L /a is not a constant, but rather depends on the frequency of the angle of 
attack variations. 

The same analysis could be applied for more complicated aerodynamic transfer functions, which allows for the 

possibility of a more general and complicated indicial function than e ~ bt . The model identification in the frequency 
domain selects the most appropriate model terms, which is the equivalent of identifying a custom indicial function 
for the unsteady aerodynamics in the frequency domain. 
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Figure 4. Unsteady aerodynamic transfer function identification 
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Figure 5. Unsteady aerodynamic transfer function model fit in the frequency domain 
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Figure 6. Bode plot of identified C L /a aerodynamic transfer function 


V. Concluding Remarks 

A new method for identifying transfer function models from noisy data was developed and demonstrated using 
simulated data. The modeling approach uses orthogonal Fourier transform modeling functions in an equation-error 
formulation in the frequency domain, along with a predicted squared error modeling metric, to identify transfer 
function model structure and estimate unknown model parameters and uncertainties. Using Fourier-based 
modulating functions extended the capability of the technique to arbitrary non-zero endpoint conditions. It was also 
shown that the technique can be applied in real time by using a recursive Fourier transform to transform time- 
domain data into the frequency domain, then applying the transfer function identification method in the same way as 
when all the data are analyzed at the same time, in post-flight batch mode. Extensions to frequency response data 
and application to state-space model identification in the frequency domain were also described. The method was 
applied to general transfer function identification and unsteady aerodynamic transfer function modeling. 

Results from simulation investigations showed that the technique described here can be used effectively for 
accurate transfer function identification. Model structure and model parameters were identified accurately at the 
same time, and model parameter uncertainties computed from the data were found to properly characterize the 
accuracy of the parameter estimates. Further work is required to evaluate the technique on real data. 

The approach described here is simple and accurate, and produces high-quality transfer function modeling results 
from noisy data. Applications include actuator modeling, “black box” transfer function modeling, unsteady 
aerodynamic modeling, and aerodynamic transfer function modeling. 
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Table 1 . Parameter estimates and uncertainties 
for a transfer function model identified from simulated data 




Orthogonal Fourier 
Transform Functions 

Orthogonal Fourier-Based 
Modulating Functions 

Model Parameter 

True Value 

Estimate ± Std. Error 

Estimate ± Std. Error 

c 0 

1.000 

1.003 ±0.020 

1.007 ±0.016 

Cl 

0.500 

0.496 ±0.008 

0.504 +0.006 

d i 

0.159 

0.159 ±0.002 

0.161 +0.002 

d 2 

0.0253 

0.0252 ±0.0003 

0.0253 ±0.0002 


Table 2. Parameter estimates and uncertainties for an 
unsteady aerodynamic transfer function model identified from simulated data 


Model Parameter 

True Value 

Estimate ± Std. Error 

Co 

2.700 

2.699 +0.011 

C 1 

2.333 

2.336 +0.009 

d l 

0.303 

0.303 ±0.002 

C L 

2.700 


a 

-5.000 


b 

3.300 
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